using Distributions
using StatsPlots
using StatsBase
using LaTeXStrings
using CSV
using DataFrames
using StatisticalRethinking
using StatisticalRethinking: link
using LinearAlgebra
using Logging
using Random
using Turing

# setting default attributes for plots
default(label=false)
Logging.disable_logging(Logging.Warn);

4.1 Why normal distributions are normal

Code 4.1

n = rand(Uniform(-1, 1), 1000, 16);
pos = sum.(eachrow(n));
density(pos)

Code 4.2

prod(1 .+ rand(Uniform(0, .1), 12))

1.9303572320698072

Code 4.3

u = Uniform(0, .1)
growth = prod.(eachrow(1 .+ rand(u, 10000, 12)));

density(growth; label="density")
# overlay normal distribution
μ = mean(growth)
σ = std(growth)
plot!(Normal(μ, σ); label="Normal")

Code 4.4

big = prod.(eachrow(1 .+ rand(Uniform(0, 0.5), 10000, 12)));
small = prod.(eachrow(1 .+ rand(Uniform(0, 0.01), 10000, 12)));
density([big, small]; layout=(2, 1))

Code 4.5

density(log.(big))

4.2 A language for describing models

Code 4.6

w = 6
n = 9
p_grid = range(0, 1; length=100)
bin_dens = [pdf(Binomial(n, p), w) for p in p_grid]
uni_dens = [pdf(Uniform(0, 1), p) for p in p_grid];
posterior = bin_dens .* uni_dens
posterior /= sum(posterior);

4.3 Gaussian model of height

Code 4.7

d = DataFrame(CSV.File("data/Howell1.csv"));

Code 4.8

describe(d)

4 rows × 7 columns

	variable	mean	min	median	max	nmissing	eltype
	Symbol	Float64	Real	Float64	Real	Int64	DataType
1	height	138.264	53.975	148.59	179.07	0	Float64
2	weight	35.6106	4.25242	40.0578	62.9926	0	Float64
3	age	29.3444	0.0	27.0	88.0	0	Float64
4	male	0.472426	0	0.0	1	0	Int64

Code 4.9

precis(d)

┌────────┬────────────────────────────────────────────────────────────┐
│  param │    mean      std     5.5%      50%    94.5%      histogram │
├────────┼────────────────────────────────────────────────────────────┤
│ height │ 138.264  27.6024  81.1086   148.59  165.735  ▁▁▁▂▂▂▂▂▂██▆▁ │
│ weight │ 35.6106  14.7192   9.3607  40.0578  54.5029  ▁▃▄▄▃▂▃▆██▅▃▁ │
│    age │ 29.3444  20.7469      1.0     27.0   66.135      █▆▆▆▆▃▃▁▁ │
│   male │  0.4724   0.4997      0.0      0.0      1.0    █▁▁▁▁▁▁▁▁▁█ │
└────────┴────────────────────────────────────────────────────────────┘

Code 4.10

d.height

544-element Vector{Float64}:
 151.765
 139.7
 136.525
 156.845
 145.415
 163.83
 149.225
 168.91
 147.955
 165.1
 154.305
 151.13
 144.78
   ⋮
 156.21
 152.4
 162.56
 114.935
  67.945
 142.875
  76.835
 145.415
 162.56
 156.21
  71.12
 158.75

Code 4.11

d2 = d[d.age .>= 18,:];

Code 4.12

plot(Normal(178, 20); xlim=(100, 250))

Code 4.13

plot(Uniform(0, 50), xlim=(-10, 60), ylim=(0, 0.03))

Code 4.14

size = 10_000
sample_μ = rand(Normal(178, 20), size)
sample_σ = rand(Uniform(0, 50), size);
prior_h = [rand(Normal(μ, σ)) for (μ, σ) in zip(sample_μ, sample_σ)];

p1 = density(sample_μ; title="μ")
p2 = density(sample_σ; title="σ")
p3 = density(prior_h; title="prior_h")

plot(p1, p2, p3, layout=@layout [a b; c])

Code 4.15

sample_μ = rand(Normal(178, 100), size)
prior_h = [rand(Normal(μ, σ)) for (μ, σ) in zip(sample_μ, sample_σ)];

density(prior_h)
vline!([0, 272])

Code 4.16

μ_list = range(150, 160; length=100)
σ_list = range(7, 9; length=100)

log_likelihood = [
    sum(logpdf(Normal(μ, σ), d2.height))
    for μ ∈ μ_list, σ ∈ σ_list
]
log_prod = log_likelihood .+ [
    logpdf(Normal(178, 20), μ) + logpdf(Uniform(0, 50), σ)
    for μ ∈ μ_list, σ ∈ σ_list
];

max_prod = maximum(log_prod)
prob = @. exp(log_prod - max_prod);

Code 4.17

# note the transposition, that's due to Julia matrix order
contour(μ_list, σ_list, prob')

Code 4.18

heatmap(μ_list, σ_list, prob')

Code 4.19

indices = collect(Iterators.product(1:length(μ_list), 1:length(σ_list)));
sample_idx = wsample(vec(indices), vec(prob), 10_000; replace=true)
sample_μ = μ_list[first.(sample_idx)]
sample_σ = σ_list[last.(sample_idx)];

Code 4.20

scatter(sample_μ, sample_σ; alpha=0.1)

Code 4.21

p1 = density(sample_μ)
p2 = density(sample_σ)
plot(p1, p2, layout=(2,1))

Code 4.22

println(hpdi(sample_μ, alpha=0.11))
println(hpdi(sample_σ, alpha=0.11))

[153.83838383838383, 155.15151515151516]
[7.262626262626263, 8.191919191919192]

Code 4.23

d3 = sample(d2.height, 20);

Code 4.24

μ_list = range(150, 170; length=200)
σ_list = range(4, 20; length=200)

log_likelihood = [
    sum(logpdf(Normal(μ, σ), d3))
    for μ ∈ μ_list, σ ∈ σ_list
]
log_prod = log_likelihood .+ [
    logpdf(Normal(178, 20), μ) + logpdf(Uniform(0, 50), σ)
    for μ ∈ μ_list, σ ∈ σ_list
]

max_prod = maximum(log_prod)
prob2 = @. exp(log_prod - max_prod)

indices = collect(Iterators.product(1:length(μ_list), 1:length(σ_list)));
sample2_idx = wsample(vec(indices), vec(prob2), 10_000; replace=true)
sample2_μ = μ_list[first.(sample2_idx)]
sample2_σ = σ_list[last.(sample2_idx)]

scatter(sample2_μ, sample2_σ; alpha=0.1)

Code 4.25

density(sample2_σ)
μ = mean(sample2_σ)
σ = std(sample2_σ)
plot!(Normal(μ, σ); label="Normal")

Code 4.26

d = DataFrame(CSV.File("data/Howell1.csv"));
d2 = d[d.age .>= 18,:];

Code 4.27

@model function model_height(height)
    μ ~ Normal(178, 20)
    σ ~ Uniform(0, 50)
    height ~ Normal(μ, σ)
end

model_height (generic function with 2 methods)

Code 4.28

m4_1 = sample(model_height(d2.height), NUTS(), 1000)

Chains MCMC chain (1000×14×1 Array{Float64, 3}):

Iterations        = 501:1:1500
Number of chains  = 1
Samples per chain = 1000
Wall duration     = 7.07 seconds
Compute duration  = 7.07 seconds
parameters        = μ, σ
internals         = lp, n_steps, is_accept, acceptance_rate, log_density, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size

Summary Statistics
  parameters       mean       std   naive_se      mcse         ess      rhat   ⋯
      Symbol    Float64   Float64    Float64   Float64     Float64   Float64   ⋯

           μ   154.6092    0.3988     0.0126    0.0113   1239.9501    0.9990   ⋯
           σ     7.7657    0.3138     0.0099    0.0092   1149.2957    0.9995   ⋯
                                                                1 column omitted

Quantiles
  parameters       2.5%      25.0%      50.0%      75.0%      97.5% 
      Symbol    Float64    Float64    Float64    Float64    Float64 

           μ   153.8444   154.3255   154.6183   154.8714   155.4424
           σ     7.1909     7.5551     7.7388     7.9676     8.4361

Code 4.29

display.(describe(m4_1; q=[0.055, 0.945]));

Summary Statistics
  parameters       mean       std   naive_se      mcse         ess      rhat   ⋯
      Symbol    Float64   Float64    Float64   Float64     Float64   Float64   ⋯

           μ   154.6092    0.3988     0.0126    0.0113   1239.9501    0.9990   ⋯
           σ     7.7657    0.3138     0.0099    0.0092   1149.2957    0.9995   ⋯
                                                                1 column omitted

Quantiles
  parameters       5.5%      94.5% 
      Symbol    Float64    Float64 

           μ   153.9801   155.2242
           σ     7.2873     8.2936

Code 4.30

init_vals = [mean(d2.height), std(d2.height)]
chain = sample(model_height(d2.height), NUTS(), 1000, init_theta=init_vals)

Chains MCMC chain (1000×14×1 Array{Float64, 3}):

Iterations        = 501:1:1500
Number of chains  = 1
Samples per chain = 1000
Wall duration     = 0.65 seconds
Compute duration  = 0.65 seconds
parameters        = μ, σ
internals         = lp, n_steps, is_accept, acceptance_rate, log_density, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size

Summary Statistics
  parameters       mean       std   naive_se      mcse        ess      rhat    ⋯
      Symbol    Float64   Float64    Float64   Float64    Float64   Float64    ⋯

           μ   154.6284    0.4095     0.0129    0.0129   813.0170    0.9999    ⋯
           σ     7.7597    0.2926     0.0093    0.0089   984.3953    0.9992    ⋯
                                                                1 column omitted

Quantiles
  parameters       2.5%      25.0%      50.0%      75.0%      97.5% 
      Symbol    Float64    Float64    Float64    Float64    Float64 

           μ   153.8513   154.3544   154.6117   154.9108   155.4373
           σ     7.2344     7.5607     7.7465     7.9384     8.3772

Code 4.31

@model function model_height(height)
    μ ~ Normal(178, 0.1)
    σ ~ Uniform(0, 50)
    height ~ Normal(μ, σ)
end

m4_2 = sample(model_height(d2.height), NUTS(), 1000)
display.(describe(m4_2; q=[0.055, 0.945]));

Summary Statistics
  parameters       mean       std   naive_se      mcse        ess      rhat    ⋯
      Symbol    Float64   Float64    Float64   Float64    Float64   Float64    ⋯

           μ   177.8637    0.1023     0.0032    0.0030   811.9349    1.0033    ⋯
           σ    24.5429    0.9562     0.0302    0.0276   834.7287    0.9991    ⋯
                                                                1 column omitted

Quantiles
  parameters       5.5%      94.5% 
      Symbol    Float64    Float64 

           μ   177.6973   178.0280
           σ    23.0719    26.0606

Code 4.32

cov(hcat(m4_1[:μ], m4_1[:σ]))

2×2 Matrix{Float64}:
 0.159072    0.00303165
 0.00303165  0.09848

Code 4.33

c = cov(hcat(m4_1[:μ], m4_1[:σ]))
cov2cor(c, diag(c))

2×2 Matrix{Float64}:
 1.0       0.193526
 0.193526  1.0

Code 4.34

# resetrange is needed due to bug in MCMCChains: https://github.com/TuringLang/MCMCChains.jl/issues/324
# once it will be fixed, direct sampling from the chain will be possible
samp_chain = sample(resetrange(m4_1), 10_000)
samp_df = DataFrame(samp_chain)
first(samp_df, 5)

5 rows × 2 columns

	μ	σ
	Float64	Float64
1	154.911	7.40381
2	155.033	7.69794
3	154.885	8.00863
4	154.094	7.39209
5	154.939	7.69093

Code 4.35

precis(samp_df)

┌───────┬─────────────────────────────────────────────────────────┐
│ param │    mean     std     5.5%     50%    94.5%     histogram │
├───────┼─────────────────────────────────────────────────────────┤
│     μ │ 154.611  0.4013  153.968  154.62  155.247  ▁▁▂▄▇▇█▇▄▂▁▁ │
│     σ │  7.7652  0.3129   7.2866  7.7433   8.2927    ▁▁▃▆█▆▅▂▁▁ │
└───────┴─────────────────────────────────────────────────────────┘

Code 4.36

data = hcat(m4_1[:μ], m4_1[:σ])
μ = mean(data, dims=1)
σ = cov(data)
mvn = MvNormal(vec(μ), σ)
post = rand(mvn, 10_000);
print(mvn)

FullNormal(
dim: 2
μ: [154.6091697298638, 7.765669889245865]
Σ: [0.15907154466725013 0.003031652659777947; 0.003031652659777947 0.09848001526298829]
)

4.4 Linear predictions

Code 4.37

d = DataFrame(CSV.File("data/Howell1.csv"));
d2 = d[d.age .>= 18,:];
# fancy way of doing scatter(d2.weight, d2.height)
@df d2 scatter(:weight, :height)

Code 4.38

Random.seed!(2971)
N = 100
a = rand(Normal(178, 20), N)
b = rand(Normal(0, 10), N);

Code 4.39

p = hline([0, 272]; ylims=(-100, 400), xlabel="weight", ylabel="hegiht")
title!(L"\beta \sim \mathcal{N}(\mu=0,\sigma=10)")

x_mean = mean(d2.weight)
xlims = extrema(d2.weight)  # getting min and max in one pass

for (α, β) ∈ zip(a, b)
    plot!(x -> α + β * (x - x_mean); xlims=xlims, c=:black, alpha=0.3)
end
display(p)

Code 4.40

b = rand(LogNormal(0, 1), 10_000)
density(b, xlims=(0, 5), bandwidth=0.1)

Code 4.41

Random.seed!(2971)
N = 100
a = rand(Normal(178, 20), N)
b = rand(LogNormal(0, 1), N);

Code 4.42

d = DataFrame(CSV.File("data/Howell1.csv"));
d2 = d[d.age .>= 18,:]
xbar = mean(d2.weight)

@model function height_regr_model(weight, height)
    a ~ Normal(178, 20)
    b ~ LogNormal(0, 1)
    μ = @. a + b * (weight - xbar)
    σ ~ Uniform(0, 50)
    height ~ MvNormal(μ, σ)
end

m4_3 = sample(height_regr_model(d2.weight, d2.height), NUTS(), 1000)
m4_3 = resetrange(m4_3);

Code 4.43

@model function height_regr_model_exp(weight, height)
    a ~ Normal(178, 20)
    log_b ~ Normal(0, 1)
    μ = @. a + exp(log_b) * (weight - xbar)
    σ ~ Uniform(0, 50)
    height ~ MvNormal(μ, σ)
end

m4_3b = sample(height_regr_model_exp(d2.weight, d2.height), NUTS(), 1000);

Code 4.44

m4_3_df = DataFrame(m4_3)
precis(m4_3_df)

┌───────┬────────────────────────────────────────────────────────┐
│ param │    mean     std     5.5%      50%    94.5%   histogram │
├───────┼────────────────────────────────────────────────────────┤
│     a │ 154.606  0.2625  154.175  154.616  155.014  ▁▁▂▅██▅▂▁▁ │
│     b │  0.9014  0.0423    0.835   0.9021   0.9665      ▁▃██▃▁ │
│     σ │  5.1071  0.1978   4.8019   5.0921   5.4553     ▁▁▅█▅▂▁ │
└───────┴────────────────────────────────────────────────────────┘

Code 4.45

round.(cov(Matrix(m4_3_df)), digits=3)

3×3 Matrix{Float64}:
  0.069  -0.0    -0.001
 -0.0     0.002  -0.0
 -0.001  -0.0     0.039

Code 4.46

p = @df d2 scatter(:weight, :height; alpha=0.3)

chain = resetrange(m4_3)
samples = sample(chain, 1000)

a_map = mean(samples[:a])
b_map = mean(samples[:b])
plot!(x -> a_map + b_map*(x-xbar))

Code 4.47

post = sample(m4_3, 1000)
post_df = DataFrame(post)
post_df[1:5,:]

5 rows × 3 columns

	a	b	σ
	Float64	Float64	Float64
1	155.241	0.939497	5.32578
2	154.435	0.890103	5.15068
3	154.393	0.95766	5.12624
4	154.644	0.975131	5.02818
5	154.719	0.833828	5.12795

Code 4.48

N = 10
dN = d2[1:N,:]

@model function height_regr_model_N(weight, height)
    a ~ Normal(178, 20)
    b ~ LogNormal(0, 1)
    m_weight = mean(weight)
    μ = @. a + b * (weight - m_weight)
    σ ~ Uniform(0, 50)
    height ~ MvNormal(μ, σ)
end

mN = sample(height_regr_model(dN.weight, dN.height), NUTS(), 1000)
mN = resetrange(mN);

Code 4.49

post = sample(mN, 20)
post_df = DataFrame(post);

xlims = extrema(d2.weight)
ylims = extrema(d2.height)
p = @df dN scatter(:weight, :height; xlims=xlims, ylims=ylims)
title!("N = $N"; xlab="weight", ylab="height")

x_mean = mean(dN.weight)
for (a, b) ∈ zip(post_df.a, post_df.b)
    plot!(x -> a + b * (x-x_mean); c="black", alpha=0.3)
end
display(p)

Code 4.50

post = sample(m4_3, 1000)
post_df = DataFrame(post)
μ_at_50 = @. post_df.a + post_df.b * (50 - xbar);

Code 4.51

density(μ_at_50; lw=2, xlab=L"\mu|weight=50")

Code 4.52

PI(μ_at_50)

2-element Vector{Float64}:
 158.5792439935818
 159.65866731906047

Code 4.53

μ = link(post_df, [:a :b], d2.weight, xbar);
μ = hcat(μ...);
Base.size(μ), μ[1:5,1]

((1000, 352), [157.56293534189655, 157.70353697789378, 156.97988672219105, 156.85965609569527, 157.2632797224988])

Code 4.54

weight_seq = 25:70
μ = link(post_df, [:a :b], weight_seq, xbar);
μ = hcat(μ...);
Base.size(μ), μ[1:5,1]

((1000, 46), [137.04072934028633, 136.64103095564798, 135.80904017330352, 136.73820023564727, 137.81265964505803])

Code 4.55

p = plot()
for i in 1:100
    scatter!(weight_seq, μ[i,:]; c=:blue, alpha=0.2)
end
display(p)

Code 4.56

μ_mean = mean.(eachcol(μ))
μ_PI = PI.(eachcol(μ))
μ_PI = vcat(μ_PI'...);

Code 4.57

@df d2 scatter(:weight, :height; alpha=0.2, xlab="weight", ylab="height")
plot!(weight_seq, [μ_mean μ_mean]; c=:black, fillrange=μ_PI, fillalpha=0.3)

Code 4.58

post = sample(m4_3, 1000)
post = DataFrame(post)

weight_seq = 25:70
μ = map(w -> post.a + post.b * (w - xbar), weight_seq)
μ = hcat(μ...)
μ_mean = mean.(eachcol(μ))
μ_CI = PI.(eachcol(μ));

Code 4.59

sim_height = simulate(post, [:a, :b, :σ], weight_seq .- xbar);
Base.size(sim_height), sim_height[1:5,1]

((1000, 46), [131.0600330744535, 131.07613382456364, 144.0054107126664, 133.0809884065339, 139.1993923582497])

Code 4.60

height_PI = PI.(eachcol(sim_height))
height_PI = vcat(height_PI'...);

Code 4.61

@df d2 scatter(:weight, :height; alpha=0.2, xlab="weight", ylab="height")
plot!(weight_seq, [μ_mean μ_mean]; c=:black, fillrange=μ_PI, fillalpha=0.3)
plot!(weight_seq, [μ_mean μ_mean]; c=:black, fillrange=height_PI, fillalpha=0.3)

Code 4.62

post = sample(m4_3, 10_000)
post = DataFrame(post)
sim_height = simulate(post, [:a, :b, :σ], weight_seq .- xbar)
height_PI = PI.(eachcol(sim_height))
height_PI = vcat(height_PI'...);

@df d2 scatter(:weight, :height; alpha=0.2, xlab="weight", ylab="height")
plot!(weight_seq, [μ_mean μ_mean]; c=:black, fillrange=μ_PI, fillalpha=0.3)
plot!(weight_seq, [μ_mean μ_mean]; c=:black, fillrange=height_PI, fillalpha=0.3)

Code 4.63

post = sample(m4_3, 1000)
post = DataFrame(post)

sim_height = [
    [
        rand(Normal(a + b * (w - xbar), σ))
        for (a, b, σ) ∈ zip(post.a, post.b, post.σ)
    ]
    for w ∈ weight_seq
]
sim_height = hcat(sim_height...)

height_PI = PI.(eachcol(sim_height));
height_PI = vcat(height_PI'...);

@df d2 scatter(:weight, :height; alpha=0.2, xlab="weight", ylab="height")
plot!(weight_seq, [μ_mean μ_mean]; c=:black, fillrange=μ_PI, fillalpha=0.3)
plot!(weight_seq, [μ_mean μ_mean]; c=:black, fillrange=height_PI, fillalpha=0.3)

4.5 Curves from lines

Code 4.64

d = DataFrame(CSV.File("data/Howell1.csv"))
scatter(d.weight, d.height; alpha=0.3)

Code 4.65

d[!, :weight_s] = standardize(ZScoreTransform, d.weight)
d[!, :weight_s2] = d.weight_s.^2;

@model function height_regr_m2(weight_s, weight_s2, height)
    a ~ Normal(178, 20)
    b1 ~ LogNormal(0, 1)
    b2 ~ Normal(0, 1)
    μ = @. a + b1 * weight_s + b2 * weight_s2
    σ ~ Uniform(0, 50)
    height ~ MvNormal(μ, σ)
end

m4_5 = sample(height_regr_m2(d.weight_s, d.weight_s2, d.height), NUTS(), 1000)
m4_5 = resetrange(m4_5)
m4_5_df = DataFrame(m4_5);

Code 4.66

precis(m4_5_df)

┌───────┬───────────────────────────────────────────────────────────┐
│ param │    mean     std     5.5%      50%    94.5%      histogram │
├───────┼───────────────────────────────────────────────────────────┤
│     a │ 146.054  0.3497  145.501  146.059  146.576  ▁▁▂▃▅██▇▄▂▁▁▁ │
│    b1 │ 21.7228  0.2941  21.2498  21.7259  22.1891    ▁▁▁▄▇█▇▄▂▁▁ │
│    b2 │ -7.8033  0.2589  -8.1943  -7.8236  -7.3754      ▁▁▂▄█▆▅▂▁ │
│     σ │  5.8137   0.179   5.5276   5.8055   6.1068   ▁▁▁▃▆██▆▄▂▁▁ │
└───────┴───────────────────────────────────────────────────────────┘

Code 4.67

Random.seed!(1)
df = sample(m4_5_df, 1000)
weight_seq = range(-2.2, 2; length=30)

mu = link(df, (r, x) -> r.a + r.b1*x + r.b2*x^2, weight_seq)
mu = hcat(mu...)
mu_mean = mean.(eachcol(mu))
mu_PI = PI.(eachcol(mu))
mu_PI = vcat(mu_PI'...)

sim_height = simulate(df, (r, x) -> Normal(r.a + r.b1*x + r.b2*x^2, r.σ), weight_seq)
sim_height = vcat(sim_height'...);

height_PI = PI.(eachcol(sim_height))
height_PI = vcat(height_PI'...);

Code 4.68

p_square = @df d scatter(:weight_s, :height; alpha=0.3, title="Square poly")
plot!(weight_seq, mu_mean; c=:black)
plot!(weight_seq, [mu_mean mu_mean]; c=:black, fillrange=mu_PI, fillalpha=0.3)
plot!(weight_seq, [mu_mean mu_mean]; c=:black, fillrange=height_PI, fillalpha=0.3)

Code 4.69

d[!, :weight_s3] = d.weight_s.^3;

@model function height_regr_m3(weight_s, weight_s2, weight_s3, height)
    a ~ Normal(178, 20)
    b1 ~ LogNormal(0, 1)
    b2 ~ Normal(0, 1)
    b3 ~ Normal(0, 1)
    μ = @. a + b1 * weight_s + b2 * weight_s2 + b3 * weight_s3
    σ ~ Uniform(0, 50)
    height ~ MvNormal(μ, σ)
end

m4_6 = sample(height_regr_m3(d.weight_s, d.weight_s2, d.weight_s3, d.height), NUTS(), 1000)
m4_6 = resetrange(m4_6)
m4_6_df = DataFrame(m4_6)
precis(m4_6_df)

┌───────┬───────────────────────────────────────────────────────────┐
│ param │    mean     std     5.5%      50%    94.5%      histogram │
├───────┼───────────────────────────────────────────────────────────┤
│     a │ 146.399  0.3158  145.888  146.405  146.898    ▁▁▁▃▅██▆▃▁▁ │
│    b1 │ 15.2346  0.4765  14.4594  15.2243  15.9857       ▁▁▁▄█▅▁▁ │
│    b2 │  -6.213  0.2561  -6.6352  -6.2127  -5.8004     ▁▁▂▅██▄▂▁▁ │
│    b3 │  3.5764  0.2239   3.2084    3.584   3.9364       ▁▂▄██▃▁▁ │
│     σ │  4.8605   0.147   4.6456   4.8546   5.0868  ▁▁▁▄▇█▇▄▂▁▁▁▁ │
└───────┴───────────────────────────────────────────────────────────┘

Random.seed!(1)
df = sample(m4_6_df, 1000)
weight_seq = range(-2.2, 2; length=30)

func = (r, x) -> r.a + r.b1*x + r.b2*x^2 + r.b3*x^3
mu = link(df, func, weight_seq)
mu = hcat(mu...)
mu_mean = mean.(eachcol(mu))
mu_PI = PI.(eachcol(mu))
mu_PI = vcat(mu_PI'...)

func = (r, x) -> Normal(r.a + r.b1*x + r.b2*x^2 + r.b3*x^3, r.σ)
sim_height = simulate(df, func, weight_seq)
sim_height = vcat(sim_height'...);

height_PI = PI.(eachcol(sim_height))
height_PI = vcat(height_PI'...);

p_cubic = @df d scatter(:weight_s, :height; alpha=0.3, title="Cubic poly")
plot!(weight_seq, mu_mean; c=:black)
plot!(weight_seq, [mu_mean mu_mean]; c=:black, fillrange=mu_PI, fillalpha=0.3)
plot!(weight_seq, [mu_mean mu_mean]; c=:black, fillrange=height_PI, fillalpha=0.3)

plot(p_square, p_cubic; layout=(1, 2))

Code 4.70 and 4.71

Looks like Julia plots don't support change of ticks proposed in the book. But much more natural way will be to remap values we're plotting back to the original scale. Example of this is below.

μ = mean(d.weight)
σ = std(d.weight)
w = @. d.weight_s * σ + μ
scatter(w, d.height; alpha=0.3)

w_s = @. weight_seq * σ + μ
plot!(w_s, mu_mean; c=:black)
plot!(w_s, [mu_mean mu_mean]; c=:black, fillrange=mu_PI, fillalpha=0.3)
plot!(w_s, [mu_mean mu_mean]; c=:black, fillrange=height_PI, fillalpha=0.3)

Code 4.72

d = DataFrame(CSV.File("data/cherry_blossoms.csv", missingstring="NA"))
precis(d)

┌────────────┬─────────────────────────────────────────────────────────┐
│      param │   mean      std    5.5%     50%    94.5%      histogram │
├────────────┼─────────────────────────────────────────────────────────┤
│       year │ 1408.0  350.885  867.77  1408.0  1948.23  ████████████▂ │
│        doy │ 104.54    6.407   94.43   105.0    115.0       ▁▂▆██▅▂▁ │
│       temp │ 6.1419   0.6636    5.15     6.1   7.2947       ▁▄▆█▄▂▁▁ │
│ temp_upper │ 7.1852   0.9929  5.8977    7.04   8.9023       ▂██▃▁▁▁▁ │
│ temp_lower │ 5.0989   0.8503  3.7876   5.145     6.37       ▁▁▁▂▇█▂▁ │
└────────────┴─────────────────────────────────────────────────────────┘

Code 4.73

d2 = d[completecases(d[!,[:doy]]),:]
d2 = disallowmissing(d2[!,[:year,:doy]])
num_knots = 15
knots_list = quantile(d2.year, range(0, 1; length=num_knots));

Code 4.74

using BSplines

basis = BSplineBasis(3, knots_list)

16-element BSplineBasis{Vector{Float64}}:
 order: 3
 breakpoints: [812.0, 1036.0, 1174.0, 1269.0, 1377.0, 1454.0, 1518.0, 1583.0, 1650.0, 1714.0, 1774.0, 1833.0, 1893.0, 1956.0, 2015.0]

Code 4.75

p1 = plot(basis)
scatter!(knots_list, repeat([1], num_knots); xlab="year", ylab="basis", legend=false)

Code 4.76

This way of calucalting bsplines is slightly slower, than shown in the book (with pre-calculated matrix) but it is much cleaner in my perspective.

You can do comparison yourself by precalculating spline matrix outside of the model and do matrix multiplication in the model instead of spline evialutaion. Example of doing this is at code block 4.79

@model function model_splines(year, doy)
    w ~ MvNormal(zeros(length(basis)), 1)
    a ~ Normal(100, 10)
    s = Spline(basis, w)
    μ = a .+ s.(year)
    σ ~ Exponential(1)
    doy ~ MvNormal(μ, σ)
end

m4_7 = sample(model_splines(d2.year, d2.doy), NUTS(0.65; init_ϵ = 9.765625e-5), 1000)

Chains MCMC chain (1000×30×1 Array{Float64, 3}):

Iterations        = 501:1:1500
Number of chains  = 1
Samples per chain = 1000
Wall duration     = 12.52 seconds
Compute duration  = 12.52 seconds
parameters        = w[1], w[2], w[3], w[4], w[5], w[6], w[7], w[8], w[9], w[10], w[11], w[12], w[13], w[14], w[15], w[16], a, σ
internals         = lp, n_steps, is_accept, acceptance_rate, log_density, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size

Summary Statistics
  parameters       mean       std   naive_se      mcse         ess      rhat   ⋯
      Symbol    Float64   Float64    Float64   Float64     Float64   Float64   ⋯

        w[1]    -0.6812    0.9172     0.0290    0.0163   1813.8789    0.9991   ⋯
        w[2]    -0.8555    0.8626     0.0273    0.0187   2705.8346    0.9990   ⋯
        w[3]     0.8252    0.8027     0.0254    0.0213   1626.8872    1.0002   ⋯
        w[4]     0.6955    0.7619     0.0241    0.0167   1808.3373    0.9990   ⋯
        w[5]     0.5630    0.7818     0.0247    0.0127   2370.8672    1.0008   ⋯
        w[6]    -1.2257    0.7915     0.0250    0.0198   2674.2308    0.9994   ⋯
        w[7]     0.2788    0.7483     0.0237    0.0189   2438.4965    0.9996   ⋯
        w[8]     1.8091    0.7712     0.0244    0.0163   2157.9901    1.0000   ⋯
        w[9]    -0.0526    0.7338     0.0232    0.0121   2400.1936    0.9991   ⋯
       w[10]     1.9439    0.8019     0.0254    0.0152   2163.3095    0.9995   ⋯
       w[11]     0.6797    0.7354     0.0233    0.0197   1598.7700    1.0000   ⋯
       w[12]     1.2764    0.7565     0.0239    0.0146   2692.3364    0.9990   ⋯
       w[13]     1.2634    0.7793     0.0246    0.0190   2460.8566    0.9994   ⋯
       w[14]    -0.7788    0.8162     0.0258    0.0157   2129.6796    0.9998   ⋯
       w[15]    -2.9017    0.8942     0.0283    0.0180   2867.7978    0.9991   ⋯
       w[16]    -2.7443    0.8905     0.0282    0.0188   2674.3391    0.9994   ⋯
           a   104.2685    0.3292     0.0104    0.0104   1152.1117    1.0032   ⋯
           σ     6.1107    0.1658     0.0052    0.0035   2781.4362    0.9991   ⋯
                                                                1 column omitted

Quantiles
  parameters       2.5%      25.0%      50.0%      75.0%      97.5% 
      Symbol    Float64    Float64    Float64    Float64    Float64 

        w[1]    -2.5031    -1.2788    -0.6836    -0.0650     1.0794
        w[2]    -2.5142    -1.4874    -0.8544    -0.2671     0.8429
        w[3]    -0.7560     0.2729     0.8494     1.3813     2.3222
        w[4]    -0.7929     0.1645     0.6935     1.1964     2.1724
        w[5]    -0.9255     0.0087     0.5661     1.0865     2.1346
        w[6]    -2.8200    -1.7612    -1.2422    -0.6835     0.3608
        w[7]    -1.2176    -0.2135     0.2979     0.7884     1.7322
        w[8]     0.2887     1.2753     1.8034     2.3243     3.3166
        w[9]    -1.4817    -0.5690    -0.0725     0.4282     1.3626
       w[10]     0.3610     1.4107     1.9837     2.5043     3.4097
       w[11]    -0.8514     0.1983     0.6660     1.1726     2.1750
       w[12]    -0.1970     0.7800     1.2584     1.8213     2.7546
       w[13]    -0.3652     0.7719     1.2357     1.8412     2.7361
       w[14]    -2.3775    -1.3282    -0.7610    -0.2397     0.8342
       w[15]    -4.7354    -3.4840    -2.8894    -2.3133    -1.1709
       w[16]    -4.5349    -3.3456    -2.7131    -2.1420    -1.0208
           a   103.6152   104.0409   104.2619   104.4908   104.8871
           σ     5.7975     5.9989     6.1047     6.2182     6.4409

Code 4.77

post = DataFrame(m4_7)

# convert columns w[*] into single column w
w_df = DataFrames.select(post, r"w")
post = DataFrames.select(post, Not(r"w"))
post[!,:w] = Vector.(eachrow(w_df))

# vector of 16 average w values
w_mean = mean.(eachcol(w_df))
p2 = plot(basis .* w_mean)
scatter!(knots_list, repeat([1], num_knots); xlab="year", ylab="basis × weight")

Code 4.78

μ = link(post, (r, x) -> r.a + Spline(basis, r.w)(x), d2.year)
μ = vcat(μ'...);

μ_PI = PI.(eachrow(μ))
μ_PI = vcat(μ_PI'...);

p3 = @df d2 scatter(:year, :doy; alpha=0.3)
μ_mean = mean.(eachrow(μ_PI))
plot!(d2.year, [μ_mean, μ_mean]; c=:black, fillrange=μ_PI, fillalpha=0.3, alpha=0)

plot(p1, p2, p3; layout=(3, 1))

Code 4.79

How to build the model with explicit spline matrix calculation

basis = BSplineBasis(3, knots_list)

# list of splines with 1 only at corresponding basis index
splines = [
    Spline(basis, [float(idx == knot) for idx ∈ 1:length(basis)])
    for knot ∈ 1:length(basis)
]

# calculate each spline for every year. Resulting matrix B is 827x16
B = [
    map(s -> s(year), splines)
    for year in d2.year
]
B = vcat(B'...);


# do not need years parameter anymore, all the information is in B matrix
@model function model_splines_matrix(doy)
    w ~ MvNormal(zeros(length(basis)), 1)
    a ~ Normal(100, 10)
    μ = a .+ B * w
    σ ~ Exponential(1)
    doy ~ MvNormal(μ, σ)
end

m4_7alt = sample(model_splines_matrix(d2.doy), NUTS(0.65; init_ϵ = 0.0001953125), 1000)

Chains MCMC chain (1000×30×1 Array{Float64, 3}):

Iterations        = 501:1:1500
Number of chains  = 1
Samples per chain = 1000
Wall duration     = 4.89 seconds
Compute duration  = 4.89 seconds
parameters        = w[1], w[2], w[3], w[4], w[5], w[6], w[7], w[8], w[9], w[10], w[11], w[12], w[13], w[14], w[15], w[16], a, σ
internals         = lp, n_steps, is_accept, acceptance_rate, log_density, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size

Summary Statistics
  parameters       mean       std   naive_se      mcse         ess      rhat   ⋯
      Symbol    Float64   Float64    Float64   Float64     Float64   Float64   ⋯

        w[1]    -0.6977    0.9024     0.0285    0.0185   2200.5259    0.9990   ⋯
        w[2]    -0.8587    0.8423     0.0266    0.0170   2112.9445    1.0009   ⋯
        w[3]     0.8351    0.7801     0.0247    0.0161   2085.0565    0.9990   ⋯
        w[4]     0.6791    0.6900     0.0218    0.0167   1563.5078    0.9990   ⋯
        w[5]     0.5912    0.8079     0.0255    0.0170   2520.6431    0.9990   ⋯
        w[6]    -1.2373    0.7507     0.0237    0.0151   2291.9338    0.9995   ⋯
        w[7]     0.2550    0.8036     0.0254    0.0215   1958.5399    0.9991   ⋯
        w[8]     1.7987    0.8424     0.0266    0.0168   1915.4226    0.9994   ⋯
        w[9]    -0.0478    0.7916     0.0250    0.0161   2057.6024    0.9993   ⋯
       w[10]     1.9124    0.7558     0.0239    0.0170   1717.8319    1.0001   ⋯
       w[11]     0.6770    0.8081     0.0256    0.0220   1902.8984    0.9991   ⋯
       w[12]     1.2699    0.7636     0.0241    0.0184   1858.3391    0.9991   ⋯
       w[13]     1.2766    0.7838     0.0248    0.0145   2455.5405    0.9990   ⋯
       w[14]    -0.7512    0.7893     0.0250    0.0141   1702.2662    0.9990   ⋯
       w[15]    -2.8963    0.8600     0.0272    0.0238   1818.3626    0.9993   ⋯
       w[16]    -2.7241    0.9117     0.0288    0.0175   2450.0885    0.9992   ⋯
           a   104.2644    0.3407     0.0108    0.0100    952.1175    1.0038   ⋯
           σ     6.0974    0.1558     0.0049    0.0038   1775.3299    0.9990   ⋯
                                                                1 column omitted

Quantiles
  parameters       2.5%      25.0%      50.0%      75.0%      97.5% 
      Symbol    Float64    Float64    Float64    Float64    Float64 

        w[1]    -2.5281    -1.2683    -0.7015    -0.1222     1.1072
        w[2]    -2.5470    -1.3879    -0.8532    -0.2840     0.8402
        w[3]    -0.6997     0.2898     0.8241     1.3810     2.3930
        w[4]    -0.6385     0.2026     0.6853     1.1511     2.0438
        w[5]    -1.0201     0.0450     0.6056     1.1317     2.2009
        w[6]    -2.6198    -1.7510    -1.2338    -0.7059     0.1801
        w[7]    -1.3066    -0.3101     0.2989     0.8075     1.7683
        w[8]     0.1302     1.2344     1.8055     2.3861     3.4671
        w[9]    -1.6704    -0.5649    -0.0236     0.4951     1.4676
       w[10]     0.4834     1.3700     1.9061     2.4198     3.3756
       w[11]    -0.8996     0.1316     0.6925     1.2162     2.2533
       w[12]    -0.1711     0.7374     1.2467     1.7692     2.7470
       w[13]    -0.1589     0.7075     1.2470     1.8217     2.8017
       w[14]    -2.2540    -1.2961    -0.7754    -0.1916     0.7812
       w[15]    -4.5643    -3.4532    -2.9137    -2.2586    -1.3194
       w[16]    -4.3992    -3.2736    -2.7774    -2.1793    -0.8404
           a   103.5640   104.0427   104.2656   104.4941   104.9266
           σ     5.8152     5.9908     6.0953     6.1964     6.4151